#ifndef _myloader_
#define _myloader_
#include "utils.h"
#include "data_readers.h"
#include "GenThreads.h"
#include <utility>
//#define VERBOSE
class CGetData{
public:
	CGetData(string file, int type):m_infile(file),
		m_type(type),pID(NULL),pHSML(NULL),  
		pRHO(NULL),  pVEL(NULL),
		pPOS(NULL),m_time(0.0)
		{    

		if(! ReadData()){
			cout<<"Fail to read data...exiting"<<endl;
			cout<<"Whereis the file?:"<<m_infile<<endl;
			EXIT_FAIL;}
		};
	~CGetData(){
		safeFree(pID);
		safeFree(pPOS);
		safeFree(pVEL);
		safeFree(pHSML);
		safeFree(pRHO);
		};
	float dist2(float *P)
		{
		return P[0]*P[0]+P[1]*P[1]+P[2]*P[2];  
		}; 
	float dist(float *P)
		{
		return sqrt(dist2(P));  
		};
	void GetCom(float *COM){
		unsigned int n=nelem,i=0;
		std::fill(COM, COM+6, 0.0f);
		for(i=0;i<n;i+=3)
			{
			COM[0]+=pPOS[i];
			COM[1]+=pPOS[i+1];
			COM[2]+=pPOS[i+2];
			//#ifdef LOAD_VEL
			COM[3]+=pVEL[i];
			COM[4]+=pVEL[i+1];
			COM[5]+=pVEL[i+2];
			//#endif
			}

		COM[0]*=1.0f/n;
		COM[1]*=1.0f/n;
		COM[2]*=1.0f/n;
		COM[3]*=1.0f/n;
		COM[4]*=1.0f/n;
		COM[5]*=1.0f/n;
		};

	void GetComByAFOF(float *COM, string afofmask);
	void GetComByMST(float *COM);
	unsigned int GetNp(){ return nelem; };
	double GetTime(){return m_time;};
	void Rotate(float phi)
		{
		};
	void Translate(const float Pos[3],const  float Vel[3])
		{
			{ 	
			//	scoped_timer timethis("Count for Transforme by loop: ");	
#pragma omp parallel for
			for(unsigned int  ic=0; ic<nelem;ic++)
				{
				unsigned int  i=3*ic;
				pPOS[i]   -= Pos[0];
				pPOS[i+1] -= Pos[1];
				pPOS[i+2] -= Pos[2];
				pVEL[i]   -= Vel[0];
				pVEL[i+1] -= Vel[1];
				pVEL[i+2] -= Vel[2];

				}  
			} 

		};
	void SetIsnap(int isnap ){m_isnap=isnap;};
protected:  
	bool ReadData()
		{
		m_gin=new CGadget(m_infile, false);
		bool ans;
		m_time=m_gin->myhead.time;
		ans = GetIDBlock(m_gin);
		if(!ans)cout<<"Problem: GetIDBlock:"<<ans<<endl;
		ans = ans && GetPosBlock(m_gin);
		if(!ans)cout<<"Problem: GetPosBlock:"<<ans<<endl;
		ans = ans && GetVelBlock(m_gin);
		if(!ans)cout<<"Problem:GetVelBlock:"<<ans<<endl;
		ans = ans && GetHsmlBlock(m_gin);
		if(!ans)cout<<"Problem:GetHsmlBlock:"<<ans<<endl;
		ans = ans && GetRhoBlock(m_gin);
		if(!ans)cout<<"Problem:GetRhoBlock:"<<ans<<endl;

		delete m_gin;
		return ans;
		};  
	bool GetRhoBlock(CGadget *g)
		{
		bool retval=true;

		nelem=0;
		safeFree(pRHO);
		nelem=g->read_block(pRHO, 
			(char *)g->rhoname[m_type].c_str(),
			m_type);

		if(nelem==0)
			{  
#ifdef VERBOSE
			cout<<"Error reading block:"<<g->rhoname[m_type].c_str()<<
				" In file: "<<g->m_filename<<endl;
#endif
			//      std::ostringstream stream;
			char buf[1024];
			sprintf(buf, "%s_rho_%d",g->m_filename.c_str(), m_type); 
			string frho=string(buf);//g->m_filename+'_rho_'+stream.str();

			//cout<<"Try to read from"<<frho<<endl;      
			CGadget *grho= new CGadget(frho, false);
			safeFree(pRHO);
			int tt=0;
			if(grho->myhead.npart[0]==0)tt=m_type;
			nelem=grho->read_block(pRHO, 
				(char *)g->rhoname[m_type].c_str(),
				tt);
			delete grho;
			retval=(bool)nelem;
			}
#ifdef VERBOSE
		cout<<"Getting: "<<g->rhoname[m_type]<<" with "<<nelem<<" elems.."<<endl;
#endif
		return retval;

		};
	bool GetHsmlBlock(CGadget *g)
		{
		bool retval=true;

		nelem=0;
		safeFree(pHSML);
		nelem=g->read_block(pHSML, 
			(char *)g->hsmlname[m_type].c_str(),
			m_type);

		if(nelem==0)
			{ 
#ifdef VERBOSE
			cout<<"Error reading block:"<<g->hsmlname[m_type].c_str()<<
				" In file: "<<g->m_filename<<endl;
			cout<<"Error reading block:"<<g->hsmlname[m_type].c_str()<<
				" In file: "<<g->m_filename<<endl;
#endif
			//      std::ostringstream stream;
			char buf[1024];
			sprintf(buf, "%s_rho_%d",g->m_filename.c_str(), m_type); 
			string frho=string(buf);//g->m_filename+'_rho_'+stream.str();
#ifdef VERBOSE
			cout<<"Try to read from"<<frho<<endl;      
#endif
			CGadget *grho= new CGadget(frho, false);
			safeFree(pHSML);
			int tt=0;
			if(grho->myhead.npart[0]==0)tt=m_type;
			nelem=grho->read_block(pHSML, 
				(char *)g->hsmlname[m_type].c_str(),
				tt);
			delete grho;
			retval=(bool)nelem;
#ifdef VERBOSE
			cout<<"Got Nelem="<<nelem<<"and retval="<<retval<<endl;
#endif
			}
#ifdef VERBOSE
		else 
			cout<<"Getting: "<<g->hsmlname[m_type]<<" with "<<nelem<<" elems.."<<endl;
#endif
		return retval;

		};
	bool GetIDBlock(CGadget *g)
		{
		bool retval=true;

		nelem=0;
		safeFree(pID);
		nelem=g->read_block(pID, 
			"ID  ",
			m_type);

		if(nelem==0)
			{  cout<<"Error reading block: "<<"ID  "<<
			" In file: "<<g->m_filename<<endl;
		retval=false;
			}
#ifdef VERBOSE
		else      cout<<"Getting: "<<"ID  "<<" with "<<nelem<<" elems.."<<endl;
#endif
		return retval;

		};


	bool GetPosBlock(CGadget *g)
		{
		bool retval=true;

		nelem=0;
		safeFree(pPOS);
		string blname = "POS ";
		nelem=g->read_blockv3(pPOS, 
			"POS ",
			m_type);

		if(nelem==0)
			{  cout<<"Error reading block:"<<"POS"<<
			" In file: "<<g->m_filename<<endl;
		retval=false;
			}
#ifdef VERBOSE
		else cout<<"Getting: "<<"POS"<<" with "<<nelem<<" elems.."<<endl;
#endif
		return retval;

		};
	bool GetVelBlock(CGadget *g)
		{
		bool retval=true;

		nelem=0;
		safeFree(pVEL);
		string blname="VEL ";
		nelem=g->read_blockv3(pVEL, 
			"VEL ",
			m_type);

		if(nelem==0)
			{  cout<<"Error reading block:"<<"VEL"<<
			" In file: "<<g->m_filename<<endl;
		retval=false;
			}
#ifdef VERBOSE
		else cout<<"Getting: "<<"VEL"<<" with "<<nelem<<" elems.."<<endl;
#endif
		return retval;

		};

protected:
	int m_type;
	double m_time;
	unsigned int nelem;
	string m_infile;
	CGadget *m_gin;
	int m_isnap;
public:
	int *pID;
	float *pPOS;
	float *pVEL;
	float *pHSML;
	float *pRHO;

	};

//int load_data(string filename, int TYPE=0);

struct CComPos
	{
	CComPos(){};
	CComPos(float Pos_[3], float Vel_[3]){
		memcpy(&Pos[0], &Pos_[0], 3*sizeof(float));
		memcpy(&Vel[0], &Vel_[0], 3*sizeof(float));
		};
	CComPos(float a,float b,float c,float d,float e,float f){
		Pos[0]=a;Pos[1]=b;Pos[2]=c;
		Vel[0]=d;Vel[1]=e;Vel[2]=f;
		};
	float Pos[3];
	float Vel[3];
	friend std::ostream& operator<<(std::ostream& os, const CComPos& e)
		{      
		os.width(16);
		os.precision(4);
		os.width(16);
		os.setf( std::ios::fixed, std::ios::floatfield ) ;
		os<<" "<<e.Pos[0];
		os.width(16);
		os<<"\t"<<e.Pos[1];
		os.width(16);
		os<<"\t"<<e.Pos[2];
		os<<std::endl;

		return os;
		}
	};
class CComVec{
public:
	CComVec(string filename):m_updated(false), m_filename(filename){
		};
	bool load_com( );
	bool   GetCom(const int isnap, CComPos &c)
		{
		if(ComVec.find(isnap)!=ComVec.end())
			{
			c=ComVec[isnap];
			return true;
			}
		m_updated=true;     
		return false;

		};
	void  PutCom(const int isnap, const CComPos c)
		{
		ComVec[isnap]=c;  
		m_updated=true;  
		};
	~CComVec()
		{

		if(m_updated)
			{
			cout<<" # Comfile was UPDATED dumping to file:  "<<m_filename<<endl;
			ofstream fo(m_filename.c_str());

			for(map<int,CComPos>::iterator it=ComVec.begin();
				it!=ComVec.end(); it++)
				{

				fo<<setprecision(4)<<std::fixed<<
					(*it).first<<setw(4)<<"  "<<
					(*it).second.Pos[0]<<setw(10)<<" "<<
					(*it).second.Pos[1]<<setw(10)<<" "<<
					(*it).second.Pos[2]<<setw(10)<<" "<<
					(*it).second.Vel[0]<<setw(10)<<" "<<
					(*it).second.Vel[1]<<setw(10)<<" "<<
					(*it).second.Vel[2]<<endl;
				}
			fo.close();

			}
		};
private:
	bool m_updated;
	string m_filename;
	map<int,CComPos> ComVec;
	};


typedef struct OnePeak{
	double pos[2];
	bool peak;
	OnePeak(double x, double y,bool pt):peak(pt){pos[0]=x;pos[1]=y;}
	OnePeak():peak(false){pos[0]=0;pos[1]=0;};
	friend std::ostream& operator<<(std::ostream& os, const OnePeak& e)
		{
		string strtype=e.peak?"Peak":"Valley";
		os<<e.pos[0]<<" "<<e.pos[1]<<" "<<
			strtype<<endl;
		return os; 
		}
	} tPeaks;

typedef vector<tPeaks> tApPeaks;
typedef vector<tPeaks>::iterator tItApPeaks;
typedef map<int,tApPeaks> tApPeaksMap;
typedef tApPeaksMap::iterator tItApPeaksMap;
class CApVec{
public:
	CApVec(string filename):m_updated(false), m_filename(filename){
		fraw=m_filename+".raw";
		load_ap();
		};
	bool load_ap();
	bool   GetAp(const int isnap, tApPeaks &c);
	void  ApCom(const int isnap, const tApPeaks c);
	void save_ap();
	~CApVec()
		{

		if(m_updated)
			{
			cout<<" # Apfile was UPDATED dumping to file:  "<<m_filename<<endl;
			save_ap();
			}
		};
	void  FillAp(int isnap);
	inline  void SetAp(vector<float> * pV){
		Ap=(*pV);  
		}
	bool AddToApPeaks(const int isnap,tApPeaks c)
		{
		bool retval=false;
		if( ApVec.find(isnap)!=ApVec.end())
			{
			retval=true;
			ApVec[isnap]=c;
			}else
				ApVec.insert(std::make_pair(isnap, c));
			save_ap();
			return retval;
		}
	vector<float> Ap;
private:
	bool m_updated;
	string m_filename;
	string fraw;
	tApPeaksMap  ApVec;

	//////////////////////////////

	};


//#undef VERBOSE
#endif

